library("here")
library("sjPlot")
library("dplyr")
library("lme4")
library("viridis")
library("lmerTest")
library("ggplot2")
library("gridExtra")
library("gt")

source("analysis_functions.R")

1. LFMM

1.1 Individual sampling

1.1.1 Summary plots

lfmm_ind <- format_lfmm(here(dirname(getwd()), "p3_methods", "outputs", "lfmm_indsampling_results.csv"))


summary_hplot(lfmm_ind, "K.1", colpal = "turbo")

summary_hplot(lfmm_ind, "TPRCOMBO", colpal = "plasma")

summary_hplot(lfmm_ind, "FDRCOMBO", colpal = "viridis", direction = -1)

1.1.2 Model summaries

run_lmer(lfmm_ind, "TPRCOMBO", table_main = "True Positive Rate")
Predictors Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp 3.35 1.12 3.00 30706 282.80 4.3 × 10−181***
sampstrat 0.46 0.15 3.00 30706 38.92 4.3 × 10−25***
K 2.54 2.54 1.00 30706 642.55 2.6 × 10−140***
m 0.62 0.62 1.00 30706 155.78 1.2 × 10−35***
phi 27.10 27.10 1.00 30706 6.86K 0.0***
H 34.89 34.89 1.00 30706 8.83K 0.0***
r 0.82 0.82 1.00 30706 206.35 1.2 × 10−46***
*** p < 0.001 
run_lmer(lfmm_ind, "FDRCOMBO", table_main = "False Discovery Rate")
Predictors Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp 1.57K 522.48 3.00 30706 3.05K 0.0***
sampstrat 33.55 11.18 3.00 30706 65.26 4.7 × 10−42***
K 15.90 15.90 1.00 30706 92.82 6.2 × 10−22***
m 29.08 29.08 1.00 30706 169.68 1.1 × 10−38***
phi 5.20 5.20 1.00 30706 30.38 3.6 × 10−8***
H 30.95 30.95 1.00 30706 180.63 4.6 × 10−41***
r 0.09 0.09 1.00 30706 0.51 0.48
*** p < 0.001

1.1.3 Megaplots

MEGAPLOT(lfmm_ind, "K.1", colpal = "turbo")

MEGAPLOT(lfmm_ind, "TPRCOMBO", colpal = "plasma")

MEGAPLOT(lfmm_ind, "FDRCOMBO", colpal = "viridis", direction = -1)

1.2 Site sampling

1.2.1 Summary plots

lfmm_site <- format_lfmm(here(dirname(getwd()), "p3_methods", "outputs", "lfmm_sitesampling_results.csv"))

summary_hplot(lfmm_site, "K.1", colpal = "turbo")

summary_hplot(lfmm_site, "TPRCOMBO", colpal = "plasma")

summary_hplot(lfmm_site, "FDRCOMBO", colpal = "viridis", direction = -1)

1.2.2 Model summaries

run_lmer(lfmm_site, "TPRCOMBO", table_main = "True Positive Rate")
Predictors Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp 0.17 0.08 2.00 17268 29.98 1.0 × 10−13***
sampstrat 0.19 0.10 2.00 17268 34.63 9.8 × 10−16***
K 0.36 0.36 1.00 17268 131.53 2.4 × 10−30***
m 0.61 0.61 1.00 17268 221.27 9.7 × 10−50***
phi 2.25 2.25 1.00 17268 818.45 6.5 × 10−176***
H 6.00 6.00 1.00 17268 2.18K 0.0***
r 0.27 0.27 1.00 17268 98.22 4.3 × 10−23***
*** p < 0.001 
run_lmer(lfmm_site, "FDRCOMBO", table_main = "False Discovery Rate")
Predictors Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp 22.32 11.16 2.00 17268 528.44 1.8 × 10−223***
sampstrat 0.43 0.22 2.00 17268 10.21 3.7 × 10−5***
K 0.23 0.23 1.00 17268 10.78 1.0 × 10−3
m 1.98 1.98 1.00 17268 93.64 4.3 × 10−22***
phi 0.50 0.50 1.00 17268 23.79 1.1 × 10−6***
H 0.36 0.36 1.00 17268 16.97 3.8 × 10−5***
r 0.11 0.11 1.00 17268 5.36 0.021**
*** p < 0.001
** p < 0.05

1.2.3 Megaplots

MEGAPLOT(lfmm_site, "K.1", colpal = "turbo")

MEGAPLOT(lfmm_site, "TPRCOMBO", colpal = "plasma")

MEGAPLOT(lfmm_site, "FDRCOMBO", colpal = "viridis", direction = -1)

2. RDA

2.1 Individual sampling

2.1.1 Summary plots

rda_ind <- format_rda(here(dirname(getwd()), "p3_methods", "outputs", "rda_indsampling_results.csv"))

summary_hplot(rda_ind, "TPRCOMBO", colpal = "plasma")

summary_hplot(rda_ind, "FDRCOMBO", colpal = "viridis", direction = -1)

2.1.2 Model summaries

run_lmer(rda_ind, "TPRCOMBO", table_main = "True Positive Rate")
Predictors Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp 0.70 0.23 3.00 30706 266.52 8.9 × 10−171***
sampstrat 0.11 0.04 3.00 30706 41.03 1.9 × 10−26***
K 0.31 0.31 1.00 30706 354.08 1.5 × 10−78***
m 0.80 0.80 1.00 30706 925.77 2.3 × 10−200***
phi 0.80 0.80 1.00 30706 919.89 4.0 × 10−199***
H 0.76 0.76 1.00 30706 879.27 1.5 × 10−190***
r 0.02 0.02 1.00 30706 21.56 3.4 × 10−6***
*** p < 0.001 
run_lmer(rda_ind, "FDRCOMBO", table_main = "False Discovery Rate")
Predictors Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp 0.04 0.01 3.00 30706 8.92 6.7 × 10−6***
sampstrat 0.00 0.00 3.00 30706 0.28 0.84
K 0.02 0.02 1.00 30706 11.89 5.6 × 10−4***
m 0.06 0.06 1.00 30706 39.03 4.2 × 10−10***
phi 0.01 0.01 1.00 30706 9.29 2.3 × 10−3**
H 0.05 0.05 1.00 30706 33.04 9.1 × 10−9***
r 0.00 0.00 1.00 30706 2.40 0.12
*** p < 0.001
** p < 0.01

1.1.3 Megaplots

MEGAPLOT(rda_ind, "TPRCOMBO", colpal = "plasma")

MEGAPLOT(rda_ind, "FDRCOMBO", colpal = "viridis", direction = -1)

2.2 Site sampling

2.2.1 Summary plots

rda_site <- format_rda(here(dirname(getwd()), "p3_methods", "outputs", "rda_sitesampling_results.csv"))

summary_hplot(rda_site, "TPRCOMBO", colpal = "plasma")

summary_hplot(rda_site, "FDRCOMBO", colpal = "viridis", direction = -1)

2.2.2 Model summaries

run_lmer(rda_site, "TPRCOMBO", table_main = "True Positive Rate")
Predictors Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp 0.13 0.07 2.00 17268 141.35 1.3 × 10−61***
sampstrat 0.02 0.01 2.00 17268 22.96 1.1 × 10−10***
K 0.07 0.07 1.00 17268 156.63 8.8 × 10−36***
m 0.13 0.13 1.00 17268 282.87 5.7 × 10−63***
phi 0.13 0.13 1.00 17268 282.87 5.7 × 10−63***
H 0.12 0.12 1.00 17268 265.42 3.1 × 10−59***
r 0.00 0.00 1.00 17268 7.18 7.4 × 10−3**
*** p < 0.001
** p < 0.01 
run_lmer(rda_site, "FDRCOMBO", table_main = "False Discovery Rate")
## boundary (singular) fit: see help('isSingular')
Predictors Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
nsamp 0.00 0.00 2.00 17270 2.00 0.140
sampstrat 0.00 0.00 2.00 17270 0.00 1.000
K 0.00 0.00 1.00 17270 0.67 0.410
m 0.00 0.00 1.00 17270 6.01 0.014***
phi 0.00 0.00 1.00 17270 6.01 0.014***
H 0.00 0.00 1.00 17270 6.01 0.014***
r 0.00 0.00 1.00 17270 0.67 0.410
*** p < 0.05

2.2.3 Megaplots

MEGAPLOT(rda_site, "TPRCOMBO", colpal = "plasma")

MEGAPLOT(rda_site, "FDRCOMBO", colpal = "viridis", direction = -1)